       PROGRAM   LSRND3D                                                00000100
CCC                                                                     00000150
CCC----LIGHT SCATTERING BY RANDOMLY ORIENTATED SPHEROIDS IN A 3-D SPACE.00000200
CCC      THE STATE OF POLARIZATION -(IL,IR,U,V)-FOR THE INCIDENT LIGHT  00000300
CCC      IS SPECIFIED IN THE PHI=0 PLANE.                               00000400
CCC      REFERENCE: S.ASANO & M. SATO (APPL.OPTICS, 1980)               00000500
CCC          NSHAPE = 1  :  PROLATE SPHEROID(IPO= 1)                    00000600
CCC          NSHAPE = 2  :  OBLATE  SPHEROID(IPO=-1)                    00000700
CCC          ZETA = INCIDENCE ANGLES                                    00000800
CCC          RFR  = COMPLEX REFRACTIVE INDX (RFRRL, RFRIM)              00000900
CCC          A = SEMI-MAJOR AXIS ; B = SEMI-MINOR AXIS ; ABYB= A/B      00001000
CCC          NMAX .LE. 50   ;   MMAX .LE. 20   ;   INCANG.LE.45         00001100
CCC          SIMPSON LAW IS USED FOR INTEGRATION OVER ANGLES            00001200
CCC          PROGRAMMED IN 1978 AT GISS(N.Y.)                           00001300
CCC-----------------------------------------MODIFIED IN 1980----------- 00001400
      IMPLICIT  REAL*8(A-H,O-Z)                                         00001500
      PARAMETER  (MQ= 70, JMX= 150)                                     00001600
      COMPLEX*16 C2,CI,CJ,RFR,RI2,DRI2,CLAMD,C2XI,DUMY1,DUMY2,          00001700
     1           DC(JMX),DI(MQ,JMX),SPCI(MQ),RE3(MQ),DRE3(MQ),RI1(MQ),  00001800
     2           DRI1(MQ),UE1(MQ,MQ),UE3(MQ,MQ),UI1(MQ,MQ),VE1(MQ,MQ),  00001900
     3           VE3(MQ,MQ),VI1(MQ,MQ),XE1(MQ,MQ),XE3(MQ,MQ),XI1(MQ,MQ),00002000
     4           YE1(MQ,MQ),YE3(MQ,MQ),YI1(MQ,MQ),F(MQ,45),G(MQ,45),    00002100
     5           ALP(MQ),BET(MQ),ALP1(MQ,21,2),BET1(MQ,21,2),           00002200
     6           DATA(2*MQ),TE1,TE2,TM1,TM2,A1(2),A2(2),A3(2),A4(2)     00002300
      COMPLEX*16 CIN(4)/(1.D0,0.D0),(0.D0,-1.D0),(-1.D0,0.D0),          00002400
     1                  (0.D0,1.D0)/                                    00002500
      COMPLEX*16 CJN(4)/(1.D0,0.D0),(0.D0,1.D0),(-1.D0,0.D0),           00002600
     2                  (0.D0,-1.D0)/                                   00002700
      DIMENSION PL(JMX),DPL(JMX),DR(JMX),DE(MQ,JMX),SIGME(MQ),SUMD(MQ), 00002800
     1          SPCE(MQ),RE1(MQ),RE2(MQ),DRE1(MQ),JWE(MQ),JWI(MQ),      00002900
     2          NUMB(21),CONV(45),RM(4,4),FD(4,4),Z(4,4,72),JWM(21,MQ), 00003000
     3          DEM(21,50,65),SIG0(MQ,45),CHI0(MQ,45),PFACT(21),        00003100
     4          EXT(45,2),SCA(45,2),P(181,4,4),SANG(181),RLTHT(181),    00003200
     5          WTHT(181),RLPHI(100),ZETA(45),DZETA(45),WZETA(45),      00003300
     6          BSCA(50,50)                                             00003400
      COMMON  /XYUV/SPCI,RI1,DRI1,RE3,DRE3,DE,RE1,DRE1,SPCE             00003500
      COMMON  /CMPLX/DI                                                 00003600
      COMMON  /BNDRY/UE1,UE3,UI1,VE1,VE3,VI1,XE1,XE3,XI1,YE1,YE3,YI1,   00003700
     1               F,G,RFR                                            00003800
      COMMON  /RMTRX/RM                                                 00003900
C                                                                       00004000
      DATA PI,CI,CJ,EPS/3.141592653589793D0,(0.D0,1.D0),(1.D0,0.D0),    00004100
     +                  1.D-10/                                         00004200
      LOGICAL*1 SHAPE                                                   00004300
      DEFINE FILE 9(1600,400,U,IV1)                                     00004400
C*     CALL  ERRSET(207,256,-1,0,0,208)                                 00004500
C                                                                       00004600
      NUNIT=9                                                           00004700
      MINT=0                                                            00004800
      MMAX=20                                                           00004900
      NLMT=50                                                           00005000
      STRD=PI/180.D0                                                    00005100
      TWOPI=2*PI                                                        00005200
      FORPI=4*PI                                                        00005300
      OMG=1.D0-EPS                                                      00005400
C*      CALL DATE(QDATE)                                                00005500
      READ(5,500) DINC,DSANG,DPHI                                       00005600
  500 FORMAT(' DINC=',F5.0,' DSANG=',F5.0,' DPHI=',F5.0)                00005700
CC          SETTING OF THE INCIDENCE ANGLES                             00005800
C           SIMPSON LAW FOR INTEGRATION OVER ANGLES                     00005900
      INCD=(90.0+EPS)/DINC                                              00006000
      INCD=2*(INCD/2)                                                   00006100
      AINCD=INCD                                                        00006200
      DINC=90.0/AINCD                                                   00006300
      INCANG=INCD+1                                                     00006400
      SIMPSN=STRD*DINC/3.                                               00006500
      WZETA(1)=SIMPSN                                                   00006600
      WZETA(INCANG)=SIMPSN                                              00006700
      DO 10 I=1,INCANG                                                  00006800
      DZETA(I)=DINC*(I-1)                                               00006900
      ZETA(I)=STRD*DZETA(I)                                             00007000
      IF(I.EQ.1.OR.I.EQ.INCANG) GO TO 15                                00007100
      WZETA(I)=2*SIMPSN                                                 00007200
      IF(MOD(I,2).EQ.0) WZETA(I)=4*SIMPSN                               00007300
   15 WZETA(I)=DSIN(ZETA(I))*WZETA(I)                                   00007400
   10 CONTINUE                                                          00007500
C          OBSERVATION POINTS (RLTHT,RLPHI)                             00007600
      NTHT0=(180.0+EPS)/DSANG                                           00007700
      NTHT0=2*(NTHT0/2)                                                 00007800
      ANTHT0=NTHT0                                                      00007900
      DSANG=180.0/ANTHT0                                                00008000
      NTHT=NTHT0+1                                                      00008100
      SIMPSN=STRD*DSANG/3.D0                                            00008200
      WTHT(1)=SIMPSN                                                    00008300
      WTHT(NTHT)=SIMPSN                                                 00008400
      DO 20 J=1,NTHT                                                    00008500
      SANG(J)=DSANG*(J-1)                                               00008600
      RLTHT(J)=STRD*SANG(J)                                             00008700
      IF(J.EQ.1.OR.J.EQ.NTHT) GO TO 25                                  00008800
      WTHT(J)=2*SIMPSN                                                  00008900
      IF(MOD(J,2).EQ.0) WTHT(J)=4*SIMPSN                                00009000
   25 WTHT(J)=DSIN(RLTHT(J))*WTHT(J)                                    00009100
   20 CONTINUE                                                          00009200
      NPHI=(360.0+EPS)/DPHI                                             00009300
      NPHI=2*(NPHI/2)                                                   00009400
      ANPHI=NPHI                                                        00009500
      NPHIS=(NPHI/2)+1                                                  00009600
      NPHI2=NPHI+2                                                      00009700
      DPHI=360.0/ANPHI                                                  00009800
      DLPHI=STRD*DPHI                                                   00009900
      DO 30 K=1,NPHI                                                    00010000
   30 RLPHI(K)=DLPHI*(K-1)                                              00010100
      DLBETA=DLPHI                                                      00010200
CC         SPHEROID  **  SIZE P.=C, A/B=ABYB, REFRACT.INDX=RFR          00010300
  400 READ(5,200,END=300) NSHAPE,ABYB,A,RFRRL,RFRIM                     00010400
  200 FORMAT(' NSHAPE=',I2,' ABYB=',F5.0,' A=',F6.0,' RFRRL=',F8.0,     00010500
     +       ' RFRIM=',F8.0 )                                           00010600
C*      WRITE(6,666) QDATE                                              00010700
C*666 FORMAT(1H1,100X,'DATE =', A8)                                     00010800
C*      CALL CLOCK                                                      00010900
      RFR=DCMPLX(RFRRL,RFRIM)                                           00011000
      SHAPE=NSHAPE.GE.2                                                 00011100
      IF(SHAPE) GO TO 11                                                00011200
C               PROLATE SPHEROIDS                                       00011300
      IPO= 1                                                            00011400
      XI=ABYB/DSQRT(ABYB*ABYB-1.D0)                                     00011500
      XI2=XI*XI-1.D0                                                    00011600
      ECCEN=1.D0/XI                                                     00011700
      C=A*ECCEN                                                         00011800
      B=C*DSQRT(XI2)                                                    00011900
      RV=(A*B*B)**(1.D0/3.D0)                                           00012000
      RTAB=DSQRT(A*A-B*B)                                               00012100 
      SAREA=0.5D0*PI*(B*B+A*A*B*DARCOS(1.D0/ABYB)/RTAB)                 00012200
      RS=DSQRT(SAREA/PI)                                                00012300
      WRITE(6,111)                                                      00012400
  111 FORMAT(1H0/20X,'*****   LIGHT SCATTERING BY PROLATE SPHEROIDAL PAR00012500
     +ICLES ORIENTED RANDOMLY IN A SPACE   *****' )                     00012600
      GO TO 22                                                          00012700
C              OBLATE SPHEROIDS                                         00012800
   11 IPO=-1                                                            00012900
      XI=1.D0/DSQRT(ABYB*ABYB-1.D0)                                     00013000
      XI2=XI*XI+1.D0                                                    00013100
      ECCEN=1.D0/DSQRT(XI2)                                             00013200
      C=A*ECCEN                                                         00013300
      B=C*XI                                                            00013400
      RV=(A*A*B)**(1.D0/3.D0)                                           00013500
      RTAB=DSQRT(A*A-B*B)                                               00013600
      SAREA=0.5D0*PI*(A*A+0.5*A*B*B/RTAB*DLOG((A+RTAB)/(A-RTAB)))       00013700
      RS=DSQRT(SAREA/PI)                                                00013800
      WRITE(6,222)                                                      00013900
  222 FORMAT(1H0/20X,'*****   LIGHT SCATTERING BY OBLATE SPHEROIDAL PART00014000
     +ICLES ORIENTED RANDOMLY IN A SPACE   *****' )                     00014100
C                                                                       00014150
   22 C1=C                                                              00014200
      C2=C1*RFR                                                         00014300
      C3=C2                                                             00014400
      WRITE(6,100) C1,RFR,ABYB,ECCEN,A,RV,B,RS                          00014500
  100 FORMAT(1H0//28X,'C   =',F8.4,16X,'REFRACTIVE INDEX =('            00014600
     + ,F8.4,',',1PE9.2,')'/28X,'A/B =',0PF7.3,17X,'ECCENTRICITY =',F8.500014700
     +  /28X,'K*A =',F8.4,16X,'EQUI-RADIUS  K*RV =',F8.5/28X,'K*B =',   00014800
     + F8.4, 16X,'EQUI-RADIUS  K*RG =',F8.5 )                           00014900
      WRITE(6,106) DINC,DSANG,DPHI                                      00015000
  106 FORMAT(1H0/28X,'DINCANG =',F6.2,'DEG',11X,'DSANG =',F6.2,'DEG',   00015100
     + 11X,'DPHI =',F6.2,'DEG'//)                                       00015200
C                                                                       00015300
CC       EMPIRICAL NUMBERS OF EXPANSION TERMS                           00015600
      IF(ABYB.GT.3.1) GO TO 40                                          00015700
      KMIN0=0.5*A+18                                                    00015800
       IF(A.LE.1.0) KMIN0= 10                                           00016000
      GO TO 50                                                          00016100
   40 KMIN0=0.7*A+30                                                    00016200
       IF(A.LT.5.0) KMIN0=0.7*A+20                                      00016300
   50 CONTINUE                                                          00016500
      MAXJ1=35+KMIN0/2+3.0*C1                                           00016600
       IF(A.LE.0.7) MAXJ1=10                                            00016700
      MAXJ2=35+KMIN0/2+3.0*CDABS(C2)                                    00016800
       IF(A.LE.0.7) MAXJ2=12                                            00016900
      IF(MAXJ1.GE.JMX) MAXJ1=JMX-2                                      00017000
      IF(MAXJ2.GE.JMX) MAXJ2=JMX-2                                      00017100
CC           INCIDENCE ANGLES -ZETA- FROM (ELV,BETA)                    00017200
      BIG=0.D0                                                          00017300
      DO 9000 I=1,INCANG                                                00017400
      DO 9010 L=1,2                                                     00017500
      EXT(I,L)=0.                                                       00017600
 9010 SCA(I,L)=0.                                                       00017700
      CONV(I)=0.D0                                                      00017800
      IF(BIG-ZETA(I)) 9020,9000,9000                                    00017900
 9020 IZETA=I                                                           00018000
      BIG=ZETA(IZETA)                                                   00018100
 9000 CONTINUE                                                          00018200
CC        AZIMUTHAL DEPENDENCE  ***  M                                  00018300
      INDEX=0                                                           00018400
      M=0                                                               00018500
 9050 CRIT0=EXT(IZETA,1)+EXT(IZETA,2)                                   00018600
      MM=M+1                                                            00018700
      M2=M+M                                                            00018800
C-    WRITE(6,102) M                                                    00018900
C-102 FORMAT(1H /11X,'*   M =',I2 )                                     00019000
      KMIN=KMIN0-DFLOAT(M)/1.3                                          00019100
C          CORRECTION FOR THE (M=0) TERM CONVERGENCE                    00019200
      IF(KMIN.LE.10) KMIN= 10                                           00019300
      IF(KMIN.GE.NLMT) KMIN=NLMT                                        00019400
      KMAX=KMIN                                                         00019500
      PFACT(MM)=1.D0                                                    00019600
      MINDX=0                                                           00019700
    8 MINDX=MINDX+1                                                     00019800
      PFACT(MM)=PFACT(MM)*(2*MINDX-1)                                   00019900
      IF(MINDX.LT.M) GO TO 8                                            00020000
CC        SPHEROIDAL RADIAL FUNCTIONS   **   SUBROUTINE                 00020100
      DO 9060 N=1,MQ                                                    00020200
      DO 9060 I=1,JMX                                                   00020300
      DE(N,I)= 0.D0                                                     00020400
 9060 DI(N,I)= (0.D0,0.D0)                                              00020500
      C1XI=-(C1*XI)**2+IPO*M*M/XI2                                      00020600
      DO 1060 NN=1,KMAX                                                 00020700
      N=M+NN-1                                                          00020800
      L=N-M                                                             00020900
      IF(SHAPE) GO TO 33                                                00021000
C                                                                       00021050
        CALL  RPRSWF(C1,XI,M,N,MINT,1,KMAX,MAXJ1,JW,DR,ALAMD,SIGME(NN), 00021100
     1             SUMD(NN),RE1(NN),DRE1(NN),RE2(NN),DRE2)              00021200
      GO TO 44                                                          00021300
   33   CALL  ROBSWF(C1,XI,M,N,MINT,1,KMAX,MAXJ1,JW,DR,ALAMD,SIGME(NN), 00021400
     1             SUMD(NN),RE1(NN),DRE1(NN),RE2(NN),DRE2)              00021500
C                                                                       00021550
   44 SPCE(NN)=ALAMD+C1XI                                               00021600
      JWE(NN)=JW                                                        00021700
      IF(MOD(L,2).EQ.0) JWE(NN)=JW+1                                    00021800
      JJA=JWE(NN)                                                       00021900
      JWM(MM,NN)=JJA                                                    00022000
      DO 1065 J=1,JJA,2                                                 00022100
      JP=(J/2)+1                                                        00022200
      DEM(MM,NN,JP)=DR(J)                                               00022300
 1065 DE(NN,J)=DR(J)                                                    00022400
      RE3(NN)=RE1(NN)+CI*RE2(NN)                                        00022500
      DRE3(NN)=DRE1(NN)+CI*DRE2                                         00022600
      NK=NN                                                             00022700
      IF(DABS(RE1(NN)).LT.1.D-70.OR.DABS(DRE2).GT.1.D+70) GO TO 1040    00022800
 1060 CONTINUE                                                          00022900
      GO TO 1050                                                        00023000
 1040 KMAX=NK                                                           00023100
      IF(KMIN.GE.KMAX) KMIN=KMAX                                        00023200
 1050 CONTINUE                                                          00023300
      C2XI=-(C2*XI)**2+IPO*M*M/XI2                                      00023400
      C3XI=-(C3*XI)**2+IPO*M*M/XI2                                      00023500
      DO 1070 NN=1,KMAX                                                 00023600
      N=M+NN-1                                                          00023700
      L=N-M                                                             00023800
      IF(RFRIM.LT.1.D-7) GO TO 1030                                     00023900
      IF(SHAPE) GO TO 55                                                00024000
C                                                                       00024050
        CALL  CPRSWF(C2,XI,M,N,MINT,2,KMAX,MAXJ2,JW,DC,CLAMD,DUMY1,     00024100
     1             DUMY2,RI1(NN),DRI1(NN),RI2,DRI2)                     00024200
       GO TO 66                                                         00024300
   55   CALL  COBSWF(C2,XI,M,N,MINT,2,KMAX,MAXJ2,JW,DC,CLAMD,DUMY1,     00024400
     1             DUMY2,RI1(NN),DRI1(NN),RI2,DRI2)                     00024500
C                                                                       00024550
   66 SPCI(NN)=CLAMD+C2XI                                               00024600
      JWI(NN)=JW                                                        00024700
      IF(MOD(L,2).EQ.0) JWI(NN)=JW+1                                    00024800
      JJB=JWI(NN)                                                       00024900
      DO 1075 J=1,JJB,2                                                 00025000
 1075 DI(NN,J)=DC(J)                                                    00025100
      GO TO 1070                                                        00025200
 1030 IF(SHAPE) GO TO 77                                                00025300
C                                                                       00025350
        CALL  RPRSWF(C3,XI,M,N,MINT,2,KMAX,MAXJ2,JW,DR,ALAMD,DUMY3,     00025400
     1              DUMY4,R1,DR1,R2,DR2)                                00025500
       GO TO 88                                                         00025600
   77   CALL  ROBSWF(C3,XI,M,N,MINT,2,KMAX,MAXJ2,JW,DR,ALAMD,DUMY3,     00025700
     1              DUMY4,R1,DR1,R2,DR2)                                00025800
C                                                                       00025850
   88 SPCI(NN)=CJ*(ALAMD+C3XI)                                          00025900
      RI1(NN)=CJ*R1                                                     00026000
      DRI1(NN)=CJ*DR1                                                   00026100
      JWI(NN)=JW                                                        00026200
      IF(MOD(L,2).EQ.0) JWI(NN)=JW+1                                    00026300
      JJB=JWI(NN)                                                       00026400
      DO 1035 J=1,JJB,2                                                 00026500
 1035 DI(NN,J)=CJ*DR(J)                                                 00026600
 1070 CONTINUE                                                          00026700
      DO 1080 NN=1,KMAX                                                 00026800
      JEN=JWE(NN)                                                       00026900
      JIN=JWI(NN)                                                       00027000
      DO 1090 LL=1,KMAX                                                 00027100
CC        INTEGRALS A-I  &  PARAMETERS U-Y                              00027200
        CALL  UVXYAI(C1,C2,XI,M,JEN,JIN,NN,LL,UE1(LL,NN),UE3(LL,NN),    00027300
     1      UI1(LL,NN),VE1(LL,NN),VE3(LL,NN),VI1(LL,NN),XE1(LL,NN),     00027400
     2      XE3(LL,NN),XI1(LL,NN),YE1(LL,NN),YE3(LL,NN),YI1(LL,NN),IPO) 00027500
CCC       PRECALCULATIONS FOR SCATTERING CROSS SECTION                  00027600
      BSCA(NN,LL)=0.D0                                                  00027700
      NAB=IABS(NN-LL)                                                   00027800
      IF(MOD(NAB,2).EQ.1) GO TO 1090                                    00027900
      DO 7020 IA=1,JEN,2                                                00028000
      I=IA-MOD(NN,2)                                                    00028100
      Q=DE(NN,IA)                                                       00028200
      IF(M.EQ.0) GO TO 7050                                             00028300
      DO 7060 IC=1,M2                                                   00028400
 7060 Q=Q*DFLOAT(I+IC)                                                  00028500
 7050 Q=2*Q*(M+I)*(M+I+1)/DFLOAT(M2+I+I+1)                              00028600
      IF(M.EQ.0) Q=2.0*Q                                                00028700
      DEAB=Q*DE(LL,IA)                                                  00028800
      IF(I.GT.NN.AND.DABS(DEAB).LT.1.D-20) GO TO 1090                   00028900
      BSCA(NN,LL)=BSCA(NN,LL)+DEAB                                      00029000
 7020 CONTINUE                                                          00029100
 1090 CONTINUE                                                          00029200
 1080 CONTINUE                                                          00029300
CC            THE OBLIQUELY INCIDENT WAVES                              00029400
      DO 3010 IN=1,INCANG                                               00029500
      ZETA0=ZETA(IN)                                                    00029600
      IF(ZETA0.LE.1.D-7) ZETA0=1.D-7                                    00029700
      ETA0=DCOS(ZETA0)                                                  00029800
      ETA20=1.0-ETA0*ETA0                                               00029900
      SQRE0=DSQRT(ETA20)                                                00030000
      PL(1)=PFACT(MM)*(SQRE0**M)                                        00030100
      PL(2)=(2*M+1)*ETA0*PL(1)                                          00030200
      DO 4030 I=1,MAXJ1                                                 00030300
      MI=M+I                                                            00030400
      FI1=I+1                                                           00030500
      PL(I+2)=((2*MI+1)*ETA0*PL(I+1)-(MI+M)*PL(I))/FI1                  00030600
 4030 DPL(I)=(MI*ETA0*PL(I)-(MI-M)*PL(I+1))/ETA20                       00030700
      DO 3040 NN=1,KMAX                                                 00030800
      N=M+NN-1                                                          00030900
      NM=N-M                                                            00031000
      JEN=JWE(NN)                                                       00031100
      NMODN=MOD(N,4)+1                                                  00031200
      FA=0.D0                                                           00031300
      GA=0.D0                                                           00031400
      S0=0.D0                                                           00031500
      DS0=0.D0                                                          00031600
      DO 3000 IR=1,JEN,2                                                00031700
      I=IR+MOD(NM,2)-1                                                  00031800
      DEPL0=DE(NN,IR)*PL(I+1)                                           00031900
      S0=S0+DEPL0                                                       00032000
      DS0=DS0+DE(NN,IR)*DPL(I+1)                                        00032100
      IF(M+I.EQ.0) GO TO 3020                                           00032200
      Q1=DE(NN,IR)/DFLOAT((M+I)*(M+I+1))                                00032300
      FA=FA+Q1*PL(I+1)                                                  00032400
      GA=GA-SQRE0*Q1*DPL(I+1)                                           00032500
      IF(I.GE.20.AND.DABS(DEPL0).LT.1.D-20) GO TO 3080                  00032600
      GO TO 3000                                                        00032700
 3020 GA= DE(NN,1)*(1.D0-ETA0)/SQRE0                                    00032800
      FA=0.D0                                                           00032900
 3000 CONTINUE                                                          00033000
 3080 Q2=4.D0                                                           00033100
      IF(M.EQ.0) Q2=2.D0                                                00033200
      F(NN,IN)=4.D0*M*CJN(NMODN)*FA/(SIGME(NN)*SQRE0)                   00033300
      G(NN,IN)=Q2*CJN(NMODN)*GA/SIGME(NN)                               00033400
      SIG0(NN,IN)=DFLOAT(M)*S0/SQRE0                                    00033500
      CHI0(NN,IN)=-SQRE0*DS0                                            00033600
 3040 CONTINUE                                                          00033700
 3010 CONTINUE                                                          00033800
CC            DETERMINATION OF EXPANSION COEFFICIENTS  -  ALP & BET     00033900
      NMAX=2*(KMIN/2)                                                   00034000
      NUMB(MM)=NMAX                                                     00034100
      DO 5000 IMODE=1,2                                                 00034200
      IF(IMODE.EQ.1) WRITE(6,103)                                       00034300
  103 FORMAT(1H ,14X,'A)  TM-MODE POLARIZATION ' )                      00034400
CC            THE BOUNDARY CONDITIONS.                                  00034500
        CALL BOUNDR(INCANG,M,KMIN,IMODE,RE1,RE2,RI1)                    00034600
      DO 3050 IN=1,INCANG                                               00034700
      DO 3030 NE=1,NMAX                                                 00034800
      N=M+NE-1                                                          00034900
      NMODN=MOD(N,4)+1                                                  00035000
      ALP(NE)= CIN(NMODN)*DI(IN,NE)                                     00035100
      BET(NE)= CIN(NMODN)*DI(IN,NE+MQ)                                  00035200
      DATA(NE)=ALP(NE)                                                  00035300
      DATA(NE+MQ)=BET(NE)                                               00035400
 3030 CONTINUE                                                          00035500
       INDEX=INDEX+1                                                    00035600
       WRITE(NUNIT'INDEX) DATA                                          00035700
CC            SCATTERING CROSS-SECTIONS                                 00035800
      IF(ZETA(IN).LT.1.D-8.AND.M.NE.1) GO TO 3050                       00035900
      IF(M.GE.3.AND.CONV(IN).LT.1.D-7) GO TO 3050                       00036000
      ASCA=0.D0                                                         00036100
      DO 7000 NA=1,NMAX                                                 00036200
      DO 7010 NB=1,NMAX                                                 00036300
      DSCA=ALP(NA)*DCONJG(ALP(NB))+BET(NA)*DCONJG(BET(NB))              00036400
      GSCA=BSCA(NA,NB)*DSCA                                             00036500
      ASCA=ASCA+GSCA                                                    00036600
 7010 CONTINUE                                                          00036700
 7000 CONTINUE                                                          00036800
      BEXT=0.D0                                                         00036900
      DO 7070 NA=1,NMAX                                                 00037000
      CEXT=ALP(NA)*SIG0(NA,IN)+BET(NA)*CHI0(NA,IN)                      00037100
      BEXT=BEXT+CEXT                                                    00037200
 7070 CONTINUE                                                          00037300
 7080 BEXT=-4.D0*BEXT                                                   00037400
      EXT(IN,IMODE)=EXT(IN,IMODE)+BEXT                                  00037500
      SCA(IN,IMODE)=SCA(IN,IMODE)+ASCA                                  00037600
      CONV(IN)=DABS(BEXT/EXT(IN,IMODE))                                 00037700
      IF(IMODE.EQ.2) GO TO 3050                                         00037800
      IF(IN.GT.1.AND.IN.LT.INCANG) GO TO 3050                           00037900
      WRITE(6,101) IN,DZETA(IN),BEXT,ASCA,NMAX                          00038000
  101 FORMAT(1H ,14X,1H(,I2,21H )   INCIDENT ANGLE =,F5.1,5H DEG. ,     00038100
     1 10X,6HAEXT =,1PD12.5,7X,6HASCA =,D12.5,5X,7H(NMAX =,I3,2H ) )    00038200
 3050 CONTINUE                                                          00038300
 5000 CONTINUE                                                          00038400
      CRIT1=EXT(IZETA,1)+EXT(IZETA,2)                                   00038500
      CRIT=DABS((CRIT1-CRIT0)/CRIT1)                                    00038600
      IF(M.GE.MMAX) GO TO 7                                             00038700
      IF(CRIT.LE.0.0001) GO TO 7                                        00038800
      M=M+1                                                             00038900
      GO TO 9050                                                        00039000
    7 MFIN=MM                                                           00039100
       REWIND NUNIT                                                     00039200
      CEXT=0.D0                                                         00039300
      CSCA=0.D0                                                         00039400
      DO 6000 IN=1,INCANG                                               00039500
      CEXT=CEXT+0.5*(EXT(IN,1)+EXT(IN,2))*WZETA(IN)                     00039600
      CSCA=CSCA+0.5*(SCA(IN,1)+SCA(IN,2))*WZETA(IN)                     00039700
 6000 CONTINUE                                                          00039800
C             AMPLITUDE FUNCTIONS                                       00039900
      DO 9070 JT=1,NTHT                                                 00040000
      DO 9080 J=1,4                                                     00040100
      DO 9080 I=1,4                                                     00040200
 9080 P(JT,I,J)=0.D0                                                    00040300
 9070 CONTINUE                                                          00040400
      DO 8080 IN=1,INCANG                                               00040500
      W=WZETA(IN)*DLBETA/FORPI                                          00040600
      CRITM=1.D0                                                        00040700
      INDEX2=IN                                                         00040800
      DO 8073 MM=1,MFIN                                                 00040900
      NMAX=NUMB(MM)                                                     00041000
      DO 8071 IMODE=1,2                                                 00041100
       READ(NUNIT'INDEX2) DATA                                          00041200
      DO 9997 NC=1,NMAX                                                 00041300
      ALP1(NC,MM,IMODE)=DATA(NC)                                        00041400
      BET1(NC,MM,IMODE)=DATA(NC+MQ)                                     00041500
 9997 CONTINUE                                                          00041600
      INDEX2=INDEX2+INCANG                                              00041700
 8071 CONTINUE                                                          00041800
 8073 CONTINUE                                                          00041900
      DO 8000 JT=1,NTHT                                                 00042000
      DO 8030 KP=1,NPHIS                                                00042100
      ISYM=2                                                            00042200
      IF(KP.EQ.1.OR.KP.EQ.NPHIS) ISYM=1                                 00042300
        CALL  ANGLE(ZETA(IN),0.D0,RLTHT(JT),RLPHI(KP),STHT,SPHI,SALPH)  00042400
      SPHI0=TWOPI-SPHI                                                  00042500
      SALPH0=TWOPI-SALPH                                                00042600
      NPHI0=NPHI2-KP                                                    00042700
       IF(JT.EQ.1. AND.KP.GE.2)  GO TO 8075                             00042800
       IF(JT.EQ.NTHT.AND.KP.GE.2)  GO TO 8075                           00042900
      DO 4050 L=1,ISYM                                                  00043000
      A1(L)=(0.D0,0.D0)                                                 00043100
      A2(L)=(0.D0,0.D0)                                                 00043200
      A3(L)=(0.D0,0.D0)                                                 00043300
 4050 A4(L)=(0.D0,0.D0)                                                 00043400
      ETA=DCOS(STHT)                                                    00043500
      ETA2=1.D0-ETA*ETA                                                 00043600
      SQRE=DSQRT(ETA2)                                                  00043700
      ABSETA=DABS(ETA)                                                  00043800
      DO 8070 MM=1,MFIN                                                 00043900
      M=MM-1                                                            00044000
      NMAX=NUMB(MM)                                                     00044100
      TE1=(0.D0,0.D0)                                                   00044200
      TE2=(0.D0,0.D0)                                                   00044300
      TM1=(0.D0,0.D0)                                                   00044400
      TM2=(0.D0,0.D0)                                                   00044500
      IF(ZETA(IN).LE.1.D-7.AND.M.NE.1) GO TO 8020                       00044600
      IF(M.GE.3.AND.CRITM.LT.1.D-6) GO TO 8020                          00044700
      IF(ABSETA.GE.OMG) GO TO 4000                                      00044800
      PL(1)=PFACT(MM)*(SQRE**M)                                         00044900
      PL(2)=(2*M+1)*ETA*PL(1)                                           00045000
      DO 4010 I=1,MAXJ1                                                 00045100
      MI=M+I                                                            00045200
      FI1=I+1                                                           00045300
      PL(I+2)=((2*MI+1)*ETA*PL(I+1)-(MI+M)*PL(I))/FI1                   00045400
 4010 DPL(I)=(MI*ETA*PL(I)-(MI-M)*PL(I+1))/ETA2                         00045500
 4000 CONTINUE                                                          00045600
      DO 8010 NC=1,NMAX                                                 00045700
      N=M+NC-1                                                          00045800
      NMMOD=MOD(N-M,2)                                                  00045900
      S=0.D0                                                            00046000
      DS=0.D0                                                           00046100
      SUM=0.D0                                                          00046200
      JEN=JWM(MM,NC)                                                    00046300
      DO 8040 IR=1,JEN,2                                                00046400
      I=IR                                                              00046500
      IF(NMMOD.EQ.0) I=IR-1                                             00046600
      IP=(I/2)+1                                                        00046700
      IF(ABSETA.GE.OMG) GO TO 8060                                      00046800
      DDPL=DEM(MM,NC,IP)*DPL(I+1)                                       00046900
      DS=DS+DDPL                                                        00047000
      IF(I.GT.NC.AND.DABS(DDPL).LT.1.D-20) GO TO 8050                   00047100
      IF(M.EQ.0) GO TO 8040                                             00047200
      S=S+DEM(MM,NC,IP)*PL(I+1)                                         00047300
      GO TO 8040                                                        00047400
 8060 SUM=SUM+.5D0*(I+1)*(I+2)*DEM(MM,NC,IP)                            00047500
 8040 CONTINUE                                                          00047600
 8050 IF(ABSETA.GE.OMG) GO TO 3                                         00047700
      SIG=M*S/SQRE                                                      00047800
      CHI=-SQRE*DS                                                      00047900
      GO TO 4                                                           00048000
    3 IF(M.NE.1) GO TO 8010                                             00048100
      IF(DABS(1.D0+ETA).LT.EPS) GO TO 6                                 00048200
      SIG=SUM                                                           00048300
      CHI=SIG                                                           00048400
      GO TO 4                                                           00048500
    6 SIG=SUM                                                           00048600
      IF(MOD(N-1,2).EQ.1) SIG=-SUM                                      00048700
      CHI=-SIG                                                          00048800
    4 TM1=TM1+ALP1(NC,MM,1)*CHI+BET1(NC,MM,1)*SIG                       00048900
      TM2=TM2+ALP1(NC,MM,1)*SIG+BET1(NC,MM,1)*CHI                       00049000
      TE1=TE1+ALP1(NC,MM,2)*SIG+BET1(NC,MM,2)*CHI                       00049100
      TE2=TE2+ALP1(NC,MM,2)*CHI+BET1(NC,MM,2)*SIG                       00049200
 8010 CONTINUE                                                          00049300
 8020 CONTINUE                                                          00049400
      DO 4060 L=1,ISYM                                                  00049500
      SPHI1=SPHI                                                        00049600
      IF(L.EQ.2) SPHI1=SPHI0                                            00049700
      FAI=M*SPHI1                                                       00049800
      SINM=DSIN(FAI)                                                    00049900
      COSM=DCOS(FAI)                                                    00050000
      IF(DABS(SINM).LT.EPS) SINM=0.D0                                   00050100
      IF(DABS(COSM).LT.EPS) COSM=0.D0                                   00050200
      A1(L)=A1(L)+TE1*COSM                                              00050300
      A2(L)=A2(L)+TM2*COSM                                              00050400
      A3(L)=A3(L)-TE2*SINM                                              00050500
 4060 A4(L)=A4(L)+TM1*SINM                                              00050600
      IF(M.EQ.0) GO TO 8070                                             00050700
      CRITM=CDABS((TM1+TM2)/(A4(1)+A2(1)))                              00050800
 8070 CONTINUE                                                          00050900
 8075 CONTINUE                                                          00051000
C         TRANSFORMATION MATRIX   --   F'                               00051100
      DO 8035 L=1,ISYM                                                  00051200
       IF(JT.EQ.1.AND.KP.GE.2) GO TO 6050                               00051300
       IF(JT.EQ.NTHT.AND.KP.GE.2) GO TO 6050                            00051400
      AR1= DREAL(A1(L))                                                 00051500
      AR2= DREAL(A2(L))                                                 00051600
      AR3= DREAL(A3(L))                                                 00051700
      AR4= DREAL(A4(L))                                                 00051800
      AI1= DIMAG(A1(L))                                                 00051900
      AI2= DIMAG(A2(L))                                                 00052000
      AI3= DIMAG(A3(L))                                                 00052100
      AI4= DIMAG(A4(L))                                                 00052200
      FD(1,1)= AR2*AR2+AI2*AI2                                          00052300
      FD(1,2)= AR3*AR3+AI3*AI3                                          00052400
      FD(1,3)= AR2*AR3+AI2*AI3                                          00052500
      FD(1,4)= AR3*AI2-AR2*AI3                                          00052600
      FD(2,1)= AR4*AR4+AI4*AI4                                          00052700
      FD(2,2)= AR1*AR1+AI1*AI1                                          00052800
      FD(2,3)= AR4*AR1+AI4*AI1                                          00052900
      FD(2,4)= AR1*AI4-AR4*AI1                                          00053000
      FD(3,1)= 2.D0*(AR2*AR4+AI2*AI4)                                   00053100
      FD(3,2)= 2.D0*(AR3*AR1+AI3*AI1)                                   00053200
      FD(3,3)= (AR2*AR1+AI2*AI1)+(AR3*AR4+AI3*AI4)                      00053300
      FD(3,4)= (AR1*AI2-AR2*AI1)+(AR3*AI4-AR4*AI3)                      00053400
      FD(4,1)= 2.D0*(AR2*AI4-AR4*AI2)                                   00053500
      FD(4,2)= 2.D0*(AR3*AI1-AR1*AI3)                                   00053600
      FD(4,3)= (AR2*AI1-AR1*AI2)+(AR3*AI4-AR4*AI3)                      00053700
      FD(4,4)= (AR2*AR1+AI2*AI1)-(AR3*AR4+AI3*AI4)                      00053800
 6050 CONTINUE                                                          00053900
C         ROTATION MATRIX   --   L(-ALPHA)                              00054000
      IF(L.EQ.2) GO TO 4070                                             00054100
      KS=KP                                                             00054200
      SALPH1=SALPH                                                      00054300
      GO TO 4080                                                        00054400
 4070 KS=NPHI0                                                          00054500
      SALPH1=SALPH0                                                     00054600
 4080   CALL  ROTAT(SALPH1)                                             00054700
      DO 5020 J=1,4                                                     00054800
      DO 5020 I=1,4                                                     00054900
      Z(I,J,KS)= 0.D0                                                   00055000
      DO 5030 K=1,4                                                     00055100
 5030 Z(I,J,KS)=Z(I,J,KS)+RM(I,K)*FD(K,J)                               00055200
 5020 CONTINUE                                                          00055300
 8035 CONTINUE                                                          00055400
 8030 CONTINUE                                                          00055500
      DO 5050 L=1,NPHI                                                  00055600
      BETA=RLPHI(L)                                                     00055700
C         ROTATION MATRIX   --   L(-BETA)                               00055800
        CALL ROTAT(BETA)                                                00055900
C         INTEGRATION OVER PARTICLE ORIENTATION                         00056000
      AZMT1=RLPHI(1)-BETA+EPS                                           00056100
      IF(AZMT1.LT.0.D0) AZMT1=AZMT1+TWOPI                               00056200
      LA=(AZMT1/DLPHI)+1                                                00056300
      IF(LA.GT.NPHI) LA=LA-NPHI                                         00056400
      AZMT2=RLPHI(1)+PI-BETA+EPS                                        00056500
      IF(AZMT2.LT.0.D0) AZMT2=AZMT2+TWOPI                               00056600
      LB=(AZMT2/DLPHI)+1                                                00056700
      IF(LB.GT.NPHI) LB=LB-NPHI                                         00056800
      DO 5070 I=1,4                                                     00056900
      IF(I.GE.3) GO TO 60                                               00057000
      JJ1=1                                                             00057100
      JJ2=2                                                             00057200
      GO TO 70                                                          00057300
   60 JJ1=3                                                             00057400
      JJ2=4                                                             00057500
   70 DO 5070 J=JJ1,JJ2                                                 00057600
      DO 5080 K=1,4                                                     00057700
 5080 P(JT,I,J)=P(JT,I,J)+W*(Z(I,K,LA)+Z(I,K,LB))*RM(K,J)               00057800
 5070 CONTINUE                                                          00057900
 5050 CONTINUE                                                          00058000
 8000 CONTINUE                                                          00058100
 8080 CONTINUE                                                          00058200
CC        AVERAGED CROSS SECTIONS FOR SCATTERING & EXTINCTION           00058300
      CSCA=PI*CSCA                                                      00058400
      CEXT=PI*CEXT                                                      00058500
      QSCA=CSCA/SAREA                                                   00058600
      QEXT=CEXT/SAREA                                                   00058700
      QABS=QEXT-QSCA                                                    00058800
      ALBEDO=QSCA/QEXT                                                  00058900
        WRITE(6,804) CEXT,CSCA,SAREA                                    00059000
  804 FORMAT(1H1/10X,'***   CROSS SECTIONS AND EFFICIENCY FACTORS AVERAG00059100
     +ED OVER RANDOM ORIENTATIONS' //16X,'CEXT =',1PD13.5,10X,'CSCA =', 00059200
     + D13.5,10X,'AREA =',D13.5)                                        00059300
        WRITE(6,805) QEXT,QSCA,QABS,ALBEDO                              00059400
  805 FORMAT(1H ,15X,'QEXT =',1PD13.5,10X,'QSCA =',D13.5,10X,'QABS =',  00059500
     + D13.5,10X,'ALBD =',D13.5 )                                       00059600
CCC       SCATTERING MATRIX FOR THE STOKES PARAMETER - (IL,IR,U,V)      00060200
       WRITE(6,800)                                                     00060300
  800 FORMAT(1H0/10X,'***   SCATTERING TRANSFORMATION MATRIX FOR THE STO00060400
     +KES PARAMETERS-(IL,IR,U,V)- FOR RANDOMLY ORIENTED SPHEROIDS' )    00060500
       WRITE(6,801)                                                     00060600
  801 FORMAT(1H0,4X,'SANG',6X,'P(1,1)',6X,'P(1,2)',6X,'P(2,1)',6X,      00060700
     1 'P(2,2)',6X,'P(3,3)',6X,'P(3,4)',6X,'P(4,3)',6X,'P(4,4)',6X,'TOTA00060800
     2L',6X,'D-POL')                                                    00060900
      DO 2040 J=1,NTHT                                                  00061000
      TOTAL=0.5D0*(P(J,1,1)+P(J,1,2)+P(J,2,1)+P(J,2,2))                 00061100
      DPOL=((P(J,2,1)+P(J,2,2))-(P(J,1,2)+P(J,1,1)))/(2.D0*TOTAL)       00061200
       WRITE(6,807) SANG(J),P(J,1,1),P(J,1,2),P(J,2,1),P(J,2,2),        00061300
     +   P(J,3,3),P(J,3,4),P(J,4,3),P(J,4,4),TOTAL,DPOL                 00061400
  807 FORMAT(1H ,F8.3,2X,1P10D12.4)                                     00061500
 2040 CONTINUE                                                          00061600
CC        LINEAR DEPOLARIZATION RATIOS                                  00061700
      WRITE(6,802)                                                      00061800
  802 FORMAT(1H0/10X,'*    THE LINEAR DEPOLARIZATION RATIOS'//5X,       00061900
     + 'SANG',7X,'DP-HV',7X,'DP-VH',7X,'DELTA' )                        00062000
      DO 2050 J=1,NTHT                                                  00062100
      DPHV=P(J,2,1)/P(J,1,1)                                            00062200
      DPVH=P(J,1,2)/P(J,2,2)                                            00062300
      DELP=2.*(P(J,1,2)+P(J,2,1))/(P(J,1,1)+P(J,1,2)+P(J,2,1)+P(J,2,2)) 00062400
        WRITE(6,806) SANG(J),DPHV,DPVH,DELP                             00062500
  806 FORMAT(1H ,2X,F8.3,1P3D12.4)                                      00062600
 2050 CONTINUE                                                          00062700
CC        NORMALIZATION OF THE SCATTERING MATRIX                        00062800
      WRITE(6,803)                                                      00062900
      WRITE(6,801)                                                      00063000
  803 FORMAT(1H0/10X,'*     NORMALIZED SCATTERING MATRIX ' )            00063100
      PNORM=0.D0                                                        00063200
      COSBAR=0.D0                                                       00063300
      DO 2020 J=1,NTHT                                                  00063400
      WW=WTHT(J)                                                        00063500
      DO 2030 M=1,4                                                     00063600
      IF(M.GE.3) GO TO 80                                               00063700
      NIN=1                                                             00063800
      NFN=2                                                             00063900
      GO TO 90                                                          00064000
   80 NIN=3                                                             00064100
      NFN=4                                                             00064200
   90 DO 2030 N=NIN,NFN                                                 00064300
 2030 P(J,M,N)=4.D0*PI*P(J,M,N)/CSCA                                    00064400
      TOTAL=0.5D0*(P(J,1,1)+P(J,1,2)+P(J,2,1)+P(J,2,2))                 00064500
      DPOL=((P(J,2,1)+P(J,2,2))-(P(J,1,2)+P(J,1,1)))/(2.D0*TOTAL)       00064600
      PNORM=PNORM+WW*TOTAL                                              00064700
      COSBAR=COSBAR+WW*DCOS(RLTHT(J))*TOTAL                             00064800
        WRITE(6,807)  SANG(J),P(J,1,1),P(J,1,2),P(J,2,1),P(J,2,2),      00064900
     +    P(J,3,3),P(J,3,4),P(J,4,3),P(J,4,4),TOTAL,DPOL                00065000
 2020 CONTINUE                                                          00065300
      COSBAR=COSBAR/PNORM                                               00065400
      PNORM=2.D0*PI*PNORM/(4.D0*PI)                                     00065500
        WRITE(6,808)  COSBAR,PNORM                                      00065600
  808 FORMAT(1H0/10X,'*     ASYMMETRY PARAMETER  -  <COS> =',1PD13.5/   00065800
     + 16X,'NORMALIZATION =',D13.5 )                                    00065900
C-      CALL CLOCK(CPTIME)                                              00066000
C-      WRITE(6,6655) CPTIME                                            00066100
 6655 FORMAT(1H0/100X,'CPU-TIME =',F10.3,'SEC' )                        00066200
      GO TO 400                                                         00066300
  300 STOP                                                              00066400
      END                                                               00066500
